home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SHOOTF.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  70 lines
  1. PROCEDURE shootf(nvar: integer; VAR v1: gln2array; VAR v2: gln1array;
  2.        delv1: gln2array; delv2: gln1array; n1,n2: integer;
  3.        x1,x2,xf,eps,h1,hmin: real; VAR f: glnvar;
  4.        VAR dv1: gln2array; VAR dv2: gln1array);
  5. (* Programs using routine SHOOTF must define the types
  6. TYPE
  7.    gln1array = ARRAY [1..n1] OF real;
  8.    gln2array = ARRAY [1..n2] OF real;
  9.    glnvar = ARRAY [1..nvar] OF real;
  10.    glnvarbynvar = ARRAY [1..nvar,1..nvar];
  11.    glindx = ARRAY [1..nvar] OF integer;
  12.    glnpbynp = glnvarbynvar;
  13. in the main routine, and set the variable kmax of ODEINT to zero. *)
  14. VAR
  15.    nok,nbad,j,iv,i: integer;
  16.    sav,det: real;
  17.    y,f1,f2: glnvar;
  18.    dfdv: glnvarbynvar;
  19.    indx: glindx;
  20. BEGIN
  21.    load1(x1,v1,y);
  22.    odeint(y,nvar,x1,xf,eps,h1,hmin,nok,nbad);
  23.    score(xf,y,f1);
  24.    load2(x2,v2,y);
  25.    odeint(y,nvar,x2,xf,eps,h1,hmin,nok,nbad);
  26.    score(xf,y,f2);
  27.    j := 0;
  28.    FOR iv := 1 TO n2 DO BEGIN
  29.       j := j+1;
  30.       sav := v1[iv];
  31.       v1[iv] := v1[iv]+delv1[iv];
  32.       load1(x1,v1,y);
  33.       odeint(y,nvar,x1,xf,eps,h1,hmin,nok,nbad);
  34.       score(xf,y,f);
  35.       FOR i := 1 TO nvar DO BEGIN
  36.          dfdv[i,j] := (f[i]-f1[i])/delv1[iv]
  37.       END;
  38.       v1[iv] := sav
  39.    END;
  40.    FOR iv := 1 TO n1 DO BEGIN
  41.       j := j+1;
  42.       sav := v2[iv];
  43.       v2[iv] := v2[iv]+delv2[iv];
  44.       load2(x2,v2,y);
  45.       odeint(y,nvar,x2,xf,eps,h1,hmin,nok,nbad);
  46.       score(xf,y,f);
  47.       FOR i := 1 TO nvar DO BEGIN
  48.          dfdv[i,j] := (f2[i]-f[i])/delv2[iv]
  49.       END;
  50.       v2[iv] := sav
  51.    END;
  52.    FOR i := 1 TO nvar DO BEGIN
  53.       f[i] := f1[i]-f2[i];
  54.       f1[i] := -f[i]
  55.    END;
  56.    ludcmp(dfdv,nvar,nvar,indx,det);
  57.    lubksb(dfdv,nvar,nvar,indx,f1);
  58.    j := 0;
  59.    FOR iv := 1 TO n2 DO BEGIN
  60.       j := j+1;
  61.       v1[iv] := v1[iv]+f1[j];
  62.       dv1[iv] := f1[j]
  63.    END;
  64.    FOR iv := 1 TO n1 DO BEGIN
  65.       j := j+1;
  66.       v2[iv] := v2[iv]+f1[j];
  67.       dv2[iv] := f1[j]
  68.    END
  69. END;
  70.